namespace Eigen {

namespace internal {

template<typename Scalar>
void
rwupdt(Matrix<Scalar, Dynamic, Dynamic>& r,
	   const Matrix<Scalar, Dynamic, 1>& w,
	   Matrix<Scalar, Dynamic, 1>& b,
	   Scalar alpha)
{
	typedef DenseIndex Index;

	const Index n = r.cols();
	eigen_assert(r.rows() >= n);
	std::vector<JacobiRotation<Scalar>> givens(n);

	/* Local variables */
	Scalar temp, rowj;

	/* Function Body */
	for (Index j = 0; j < n; ++j) {
		rowj = w[j];

		/* apply the previous transformations to */
		/* r(i,j), i=0,1,...,j-1, and to w(j). */
		for (Index i = 0; i < j; ++i) {
			temp = givens[i].c() * r(i, j) + givens[i].s() * rowj;
			rowj = -givens[i].s() * r(i, j) + givens[i].c() * rowj;
			r(i, j) = temp;
		}

		/* determine a givens rotation which eliminates w(j). */
		givens[j].makeGivens(-r(j, j), rowj);

		if (rowj == 0.)
			continue; // givens[j] is identity

		/* apply the current transformation to r(j,j), b(j), and alpha. */
		r(j, j) = givens[j].c() * r(j, j) + givens[j].s() * rowj;
		temp = givens[j].c() * b[j] + givens[j].s() * alpha;
		alpha = -givens[j].s() * b[j] + givens[j].c() * alpha;
		b[j] = temp;
	}
}

} // end namespace internal

} // end namespace Eigen
